#---------------------------------#
# Conflict analysis between state #
# and groups                      #
#---------------------------------#
rm(list = ls())

war.data <- readRDS("data/war_data_replication.rds")
load("data/mids.RData")
load("data/fmids.RData")
load("data/cinc_ratio.RData")
load("data/contig.RData")
load("data/major_power.RData")
load("data/joint_democ.RData")
load("data/joint_autoc.RData")
load("data/out_group.RData")
intl_groups <- readRDS("intl_order/international_order_replication.rds")
coop_regime <- dplyr::select(
  intl_groups, 
  c(
    ccode, year, interlayer
  )
)


set.seed(1)
war.data$fatal_mids_num <- war.data$fatal_mid
war.data$fatal_mids_num <- as.numeric(as.character(war.data$fatal_mids_num))
war.data$year_num <- war.data$year - mean(war.data$year)
war.data$state_count <- as.numeric(war.data$state_count)
training_rows <- sample(nrow(war.data), round(0.75 * nrow(war.data)))
war.data.training <- war.data[training_rows, ]
war.data.testing  <- war.data[-training_rows, ]

war.data.training <- filter(
  war.data.training, 
  !is.na(fatal_mid)
)

war.data.testing <- filter(
  war.data.training, 
  !is.na(fatal_mid)
)

bivariate_regression <- fixest::feglm(fatal_mid ~ out_group | ccode1 + ccode2 + dyad + year, data = war.data.training, family = "binomial")
summary(bivariate_regression, cluster = c("year", "dyad", "ccode1", "ccode2"))

bivariate_lpm <- fixest::feols(fatal_mid ~ out_group | ccode1 + ccode2  + dyad + year, data = war.data.training)
summary(bivariate_lpm, cluster = c("year", "dyad", "ccode1", "ccode2"))

full_regression <- fixest::feglm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio | ccode1 + ccode2  + dyad + year, data = war.data.training, family = "binomial")
full_lpm <- fixest::feols(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio | ccode1 + ccode2  + dyad + year, data = war.data.training)
summary(full_regression, cluster = c("year", "dyad", "ccode1", "ccode2"))
summary(full_lpm, cluster = c("year", "dyad", "ccode1", "ccode2"))


pred_model_fe <- data.frame(
  out_group = c(0, 1, 0, 0),
  democ_joint = c(0, 0, 0, 1),
  autoc_joint = median(war.data.training$autoc_joint, na.rm = T), 
  cinc_ratio = median(war.data.training$cinc_ratio, na.rm = T),
  major_power_joint = median(war.data.training$major_power_joint, na.rm = T),
  ccode1 = "2",
  ccode2 = "740",
  dyad = "2-740",
  year = "1941"
)

check_data <- predict(full_regression, newdata = pred_model_fe, type = "response")

round(check_data[2] - check_data[1], 2)
round(check_data[4] - check_data[3], 2)

pred_model_fe <- data.frame(
  out_group = median(war.data.training$out_group, na.rm = T),
  democ_joint = median(war.data.training$democ_joint, na.rm = T),
  autoc_joint = median(war.data.training$autoc_joint, na.rm = T), 
  cinc_ratio = c(0, 1),
  major_power_joint = median(war.data.training$major_power_joint, na.rm = T),
  ccode1 = "2",
  ccode2 = "740",
  dyad = "2-740",
  year = "1941"
)

check_data <- predict(full_regression, newdata = pred_model_fe, type = "response")

round(check_data[2] - check_data[1], 2)

summary(full_regression, cluster = c("year", "dyad", "ccode1", "ccode2"))
summary(full_lpm, cluster = c("year", "dyad", "ccode1", "ccode2"))

save(full_regression, file = "model_output/full_regression.rda")
save(bivariate_regression, file = "model_output/simple_regression.rda")
save(full_lpm, file = "model_output/full_lpm.rda")

full_regression_mid <- fixest::feglm(mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio | ccode1 + ccode2  + dyad + year, data = war.data, family = "binomial")
full_lpm_mid        <- fixest::feols(mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio | ccode1 + ccode2  + dyad + year, data = war.data)

full_glm        <- glm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + contig + poly(year, 2), data = war.data, family = "binomial")
full_lm         <- lm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + contig + poly(year, 2), data = war.data)

full_glm_mid    <- glm(mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + contig + poly(year, 2), data = war.data, family = "binomial")
full_lm_mid     <- lm(mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + contig + poly(year, 2), data = war.data)

summary(full_regression_mid, cluster = c("year", "ccode1", "ccode2", "nameb"))
summary(full_lpm_mid, cluster = c("year", "dyad", "ccode1", "ccode2"))
summary(full_glm)
summary(full_lm)
summary(full_glm_mid)
summary(full_lm_mid)


save(full_regression_mid, file = "model_output/full_regression_mid.rda")
save(full_lpm_mid, file = "model_output/simple_regression_mid.rda")
save(full_glm, file = "model_output/full_glm.rda")
save(full_lm, file = "model_output/full_lm.rda")
save(full_glm_mid, file = "model_output/full_glm_mid.rda")
save(full_lm_mid, file = "model_output/full_lm_mid.rda")

reduced_regression <- fixest::feglm(fatal_mid ~ democ_joint + autoc_joint + major_power_joint + cinc_ratio | ccode1 + ccode2  + dyad + year, data = war.data.training, family = "binomial")
pred_auc_full <- c()
pred_auc_reduced <- c()
set.seed(1)
for (i in 1:5000)
{
  cat("\r")
  cat(paste0("\rPct of boot AUC: ", round(i/5000, 2)))
  boot_test <- war.data.testing[sample(nrow(war.data.testing), nrow(war.data.testing), replace = T),]
  suppressMessages({
    full_fit_full <- predict(full_regression, newdata = boot_test, type = "response")
    full_fit_redu <- predict(reduced_regression, newdata = boot_test, type = "response")
    pred_auc_full[i]    <- as.numeric(pROC::auc(boot_test$fatal_mid, full_fit_full))
    pred_auc_reduced[i] <- as.numeric(pROC::auc(boot_test$fatal_mid, full_fit_redu))
  })
}

boot_auc_conf <- data.frame(
  pred_auc_full,
  pred_auc_reduced
)

t.test.auc.conf <- t.test(pred_auc_full, pred_auc_reduced)
saveRDS(boot_auc_conf, file = "model_output/boot_auc_values_conf.rds")
saveRDS(t.test.auc.conf, file = "model_output/t.test.auc.conf.rds")


fit_full <- predict(full_regression, newdata = war.data.testing, type = "response")
fit_redu <- predict(reduced_regression, newdata = war.data.testing, type = "response")


table(round(fit_full), war.data.testing$fatal_mid)
table(round(fit_redu), war.data.testing$fatal_mid)

round(sum(diag(prop.table(table(round(fit_full), war.data.testing$fatal_mid)))), 2)
round(sum(diag(prop.table(table(round(fit_redu), war.data.testing$fatal_mid)))), 2)

tp_full <- table(round(fit_full), war.data.testing$fatal_mid)[2,2]
tp_redu <- table(round(fit_redu), war.data.testing$fatal_mid)[2,2]
fp_full <- table(round(fit_full), war.data.testing$fatal_mid)[2,1]
fp_redu <- table(round(fit_redu), war.data.testing$fatal_mid)[2,1]
fn_full <- table(round(fit_full), war.data.testing$fatal_mid)[1,2]
fn_redu <- table(round(fit_redu), war.data.testing$fatal_mid)[1,2]
precision_full <- tp_full / (tp_full + fp_full)
precision_redu <- tp_redu / (tp_redu + fp_redu)
recall_full    <- tp_full / (tp_full + fn_full)
recall_redu    <- tp_redu / (tp_redu + fn_redu)
f1_score_full <- round((2 * precision_full * recall_full)/(precision_full + recall_full), 2)
f1_score_redu <- round((2 * precision_redu * recall_redu)/(precision_redu + recall_redu), 2)

set.seed(1)
lambda_grid <- seq(0, 1, 0.25)
alpha_grid <- seq(1)
srchGrid <- expand.grid(.alpha = alpha_grid, .lambda = lambda_grid)
war.data.enet <- dplyr::select(
  war.data,
  c(
    fatal_mid, out_group, democ_joint, autoc_joint, major_power_joint, cinc_ratio, contig, year
  )
) 
war.data.enet <- war.data.enet[complete.cases(war.data.enet),]
war.data.enet$fatal_mid <- factor(war.data.enet$fatal_mid)
conf_elnet = train(
  fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + contig + poly(year, 2), 
  data = war.data.enet,
  method = "glmnet",
  tuneGrid = srchGrid,
  allowParallel = TRUE
)
save(conf_elnet, file = "model_output/conf_elnet.rda")


## TERGM analysis 
democ_joint        <- joint_democ
autoc_joint        <- joint_autoc
major_power_joint  <- major_power
fmids_out_sample   <- fmids[195:length(fmids)]
fmids              <- fmids[1:194]
out_group          <- out_group[1:194]
major_power_joint  <- major_power_joint[1:194]
contig             <- contig[1:194]
democ_joint        <- democ_joint[1:194]
autoc_joint        <- autoc_joint[1:194]
cinc_ratio         <- cinc_ratio[1:194]
# 
cores <- parallel::detectCores() - 2
cl <- parallel::makeCluster(cores, type = "PSOCK")
set.seed(1)
suppressMessages(
  btergm_output <- btergm(
    fmids ~ edges + 
      timecov(transform = function(t) t + t^2) + 
      kstar(2:3) + 
      triangle +   
      edgecov(out_group) + 
      edgecov(major_power_joint) + 
      edgecov(contig)      +
      edgecov(democ_joint) + 
      edgecov(autoc_joint) + 
      edgecov(cinc_ratio) , 
    R = 500,
    parallel = "snow", 
    ncpus = cores, 
    cl = cl
  )
)

parallel::stopCluster(cl)

summary(btergm_output, level = 0.95)
save(btergm_output, file = "model_output/btergm_output.RData")      
load("model_output/btergm_output.RData")      


cl <- makeCluster(cores, type = "PSOCK")
set.seed(1)
gof_btergm <- gof(
  btergm_output,
  target = fmids_out_sample,
  nsim = 100,
  MCMC.interval = 1000,
  MCMC.burnin = 10000,
  ncpus = cores,
  cl = cores,
  statistics = c(deg, triad.undirected, kstar)
)
parallel::stopCluster(cl)

png("data_visualization/fatal_mids_plot.png", width = 600, height = 300)
plot(gof_btergm)
dev.off()

load("data/cinc_ratio.RData")
load("data/contig.RData")
load("data/major_power.RData")
load("data/joint_democ.RData")
load("data/joint_autoc.RData")
load("data/out_group.RData")

democ_joint        <- joint_democ
autoc_joint        <- joint_autoc
major_power_joint  <- major_power
cores <- parallel::detectCores() - 2
cl <- makeCluster(cores, type = "PSOCK")
set.seed(123)
suppressMessages(
  btergm_output_mids <- btergm(
    mids ~ edges + 
      timecov(transform = function(t) t + t^2) + 
      kstar(2:3) + 
      triangle +   
      edgecov(out_group) + 
      edgecov(major_power_joint) + 
      edgecov(contig)      +
      edgecov(democ_joint) + 
      edgecov(autoc_joint) + 
      edgecov(cinc_ratio), 
    R = 500,
    parallel = "snow", 
    ncpus = cores, 
    cl = cl
  )
)

stopCluster(cl)
# 
summary(btergm_output_mids, level = 0.95)
save(btergm_output_mids, file = "model_output/btergm_output_mids.RData")      


pred_model_fe <- bind_rows(
  data.frame(
    out_group = c(0, 1),
    democ_joint = mean(as.numeric(as.character(war.data.training$democ_joint, na.rm = T))),
    autoc_joint = mean(as.numeric(as.character(war.data.training$autoc_joint, na.rm = T))),
    cinc_ratio = mean(as.numeric(as.character(war.data.training$cinc_ratio, na.rm = T))),
    major_power_joint = mean(as.numeric(as.character(war.data.training$major_power_joint, na.rm = T))),
    ccode1 = "2",
    ccode2 = "740",
    dyad = "2-740",
    year = "1941",
    label = c("In-cooperative\nregime relation", "Out-cooperative\nregime relation")
  ), 
  data.frame(
    democ_joint = c(0, 1),
    out_group =  mean(as.numeric(as.character(war.data.training$out_group, na.rm = T))),
    autoc_joint = mean(as.numeric(as.character(war.data.training$autoc_joint, na.rm = T))),
    cinc_ratio = mean(as.numeric(as.character(war.data.training$cinc_ratio, na.rm = T))),
    major_power_joint = mean(as.numeric(as.character(war.data.training$major_power_joint, na.rm = T))),
    ccode1 = "2",
    ccode2 = "740",
    dyad = "2-740",
    year = "1941",
    label = c("Not joint\ndemocracies", "Joint\ndemocracies")
  )
)

set.seed(1)
pred_values <- data.frame()
boots <- 1000
for (i in 1:boots)
{
  sample_rows <- sample(nrow(war.data.training), nrow(war.data.training), replace = T)
  
  full_data <- data.frame(
    row = sample_rows,
    weight = 1
  ) %>% aggregate(
    weight ~ row, data = ., FUN = "sum"
  ) 
  
  all_rows <- data.frame(
    row = 1:nrow(war.data.training)
  ) %>% left_join(
    full_data,
    by = "row"
  )
  all_rows$weight[is.na(all_rows$weight)] <- 0
  
  suppressWarnings({suppressMessages({boot_mod <- fixest::feglm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio | ccode1 + ccode2  + dyad + year, data = war.data.training, family = "binomial", weights = all_rows$weight)})})
  
  pred_values <- bind_cols(
    pred_model_fe, 
    data.frame(
      pred = suppressWarnings({suppressMessages({predict(boot_mod, newdata = pred_model_fe)})}),
      boot = i
    )
  ) %>% bind_rows(
    pred_values
  )
  
  cat(paste0("\rPct: ", round(i/boots, 2) * 100, "%"))
}

saveRDS(pred_values, file = "model_output/pred_values_outgroup.rds")

## Latent factor model with random year, dyad and node
war.data.training$out_group          <- as.factor(as.character(war.data.training$out_group))
war.data.training$major_power_joint  <- as.factor(as.character(war.data.training$major_power_joint))
war.data.training$contig             <- as.factor(as.character(war.data.training$contig))
war.data.training$democ_joint        <- as.factor(as.character(war.data.training$democ_joint))
war.data.training$autoc_joint        <- as.factor(as.character(war.data.training$autoc_joint))
war.data.training$year               <- as.factor(as.character(war.data.training$year))
war.data.training$namea              <- as.factor(as.character(war.data.training$ccode1))
war.data.training$nameb              <- as.factor(as.character(war.data.training$ccode2))
war.data.training$dyad               <- as.factor(as.character(war.data.training$dyad))
war.data.training$state_count        <- as.factor(as.character(war.data.training$state_count))
# 
set.seed(1)
system.time(multi_factor_model <- brm(data = war.data.training, 
                                     family = bernoulli,
                                     fatal_mid ~ 1 +
                                       (1|year) + (1|state_count) + 
                                       (1|dyad) + (1|namea) + (1|nameb) + 
                                       out_group + major_power_joint + contig + democ_joint + 
                                       autoc_joint + cinc_ratio + poly(year_num, 2), 
                                     prior = c(set_prior("normal(0, 1)", class = "b"),
                                               set_prior("cauchy(0, 1)", class = "sd")),
                                     chains = 8, iter = 2000, warmup = 1500, cores = 8,
                                     control = list(adapt_delta = 0.95)))
summary(multi_factor_model)
save(multi_factor_model, file = "model_output/multi_factor_model.RData")
multi_sum <- summary(multi_factor_model)
save(multi_sum, file = "model_output/multi_factor_summary.RData")
load("model_output/multi_factor_model.RData")
summary(multi_factor_model)


war.data.testing$out_group          <- as.factor(as.character(war.data.testing$out_group))
war.data.testing$major_power_joint  <- as.factor(as.character(war.data.testing$major_power_joint))
war.data.testing$contig             <- as.factor(as.character(war.data.testing$contig))
war.data.testing$democ_joint        <- as.factor(as.character(war.data.testing$democ_joint))
war.data.testing$autoc_joint        <- as.factor(as.character(war.data.testing$autoc_joint))
war.data.testing$year               <- as.factor(as.character(war.data.testing$year))
war.data.testing$namea              <- as.factor(as.character(war.data.testing$ccode1))
war.data.testing$nameb              <- as.factor(as.character(war.data.testing$ccode2))
war.data.testing$dyad               <- as.factor(as.character(war.data.testing$dyad))
#war.data.testing$lagged_fatal_mid   <- as.factor(as.character(war.data.testing$lagged_fatal_mid))
war.data.testing$state_count        <- as.factor(as.character(war.data.testing$state_count))

predicted_lf_list <- list()
range.1 <- seq(1, nrow(war.data.testing), 10000)
range.2 <- c(seq(10000, nrow(war.data.testing), 10000), nrow(war.data.testing))

range <- 10000
range.1 <- seq(1, nrow(war.data.testing), range)
range.2   <- unique(c(seq(1 + range - 1, nrow(war.data.testing), range), nrow(war.data.testing)))

for (i in 1:length(range.1))
{ 
 predicted_lf <- predict(multi_factor_model, type = "response", war.data.testing[range.1[i]:range.2[i],], allow_new_levels = T)
 predicted_lf_list[[i]] <- data.frame(
   predicted_lf = predicted_lf[,1],
   lf_result = war.data.testing[range.1[i]:range.2[i],]$fatal_mid
 )  
} 
save(predicted_lf_list, file = "model_output/predicted_lf_list.rda")
load("model_output/predicted_lf_list.rda")

load("data/mids.RData")
load("data/fmids.RData")
load("data/cinc_ratio.RData")
load("data/contig.RData")
load("data/major_power.RData")
load("data/joint_democ.RData")
load("data/joint_autoc.RData")
load("data/out_group.RData")



## TERGM analysis 
democ_joint        <- joint_democ
autoc_joint        <- joint_autoc
major_power_joint  <- major_power

require("parallel")
load("model_output/btergm_output.RData") 

cores <- detectCores() - 1
cl <- makeCluster(cores, type = "PSOCK")
set.seed(123)
btergm_test <- btergm(
  fmids ~ edges + 
    timecov(transform = function(t) t + t^2) + 
    kstar(2:3) + 
    triangle +     
    edgecov(out_group) + 
    edgecov(major_power_joint) + 
    edgecov(contig)      +
    edgecov(democ_joint) + 
    edgecov(autoc_joint) + 
    edgecov(cinc_ratio), 
  R = 500,
  parallel = "snow", 
  ncpus = cores, 
  cl = cl
)
summary(btergm_test)
stopCluster(cl)

btergm_output@coef
data_tergm <- btergm_test@effects
data_tergm$y <- btergm_test@response
data_tergm <- filter(
  data_tergm, 
  `edgecov.timecov1[[i]]` >= 185^2
)

inverse.logit <- function(x){
  if(!is.numeric(x)){return("Input is not numeric.")}
  y <- (1+exp(-x))^-1
  return(y)
}

y <- data_tergm$y
data_tergm <- dplyr::select(data_tergm, -y)
full_data_tergm <-  as.matrix(data_tergm)
predict_tergm <- inverse.logit(full_data_tergm %*% as.vector(btergm_output@coef))
predicted_tergm <- data.frame(
  predict_tergm = predict_tergm, 
  tergm_result = y
)

predicted_tergm <- dplyr::rename(
  predicted_tergm, 
  pred = predict_tergm, 
  result = tergm_result
) %>% 
  mutate(
    model = "TERGM"
  )

predicted_mlf <- do.call(
 bind_rows, predicted_lf_list
) %>% 
 dplyr::rename(
   pred =  predicted_lf,
   result = lf_result
 ) %>% 
 mutate(
   model = "Multiplicative\nlatent factor model"
 )

predicted_fe <- data.frame(
  pred = predict(full_regression, newdata = war.data.testing),
  result = war.data.testing$fatal_mid,
  model = "Fixed effect"
)

predicted_lpm <- data.frame(
  pred = predict(full_lpm, newdata = war.data.testing),
  result = war.data.testing$fatal_mid,
  model = "LPM"
)

all_fits <- bind_rows(
  predicted_fe, 
  predicted_tergm,
  predicted_mlf,
  predicted_lpm
) %>% 
  mutate(
    dv = "Fatal MIDs"
  )

save(all_fits, file = "model_output/conflict_model_fits.Rda")

dyad_data <- create_dyadyears(directed = F) %>% 
  add_nmc() %>% 
  add_peace_years() %>% 
  add_cow_mids() %>% 
  add_cow_majors() %>%
  add_democracy() %>% 
  add_fpsim() %>%
  transform(
    cinc_ratio = if_else(cinc1 > cinc2, (cinc1 - cinc2)/(cinc1 + cinc2), (cinc2 - cinc1)/(cinc1 + cinc2)),
    democ_joint = if_else(polity21 >= 6 & polity22 >= 6, 1, 0),
    autoc_joint = if_else(polity21 <= -6 & polity22 <= -6, 1, 0),
    dyad = paste0(ccode1, "-", ccode2),
    mid = cowmidongoing, 
    fatal_mid = if_else(cowmidongoing == 1 & fatality > 0, 1, 0),
    major_power_joint = if_else(cowmaj1 == 1 & cowmaj2 == 1, 1, 0)
  ) %>% left_join(
    coop_regime, 
    by = c(
      "ccode1" = "ccode",
      "year"
    )
  ) %>% left_join(
    coop_regime, 
    by = c(
      "ccode2" = "ccode",
      "year"
    )
  ) %>% mutate(
    out_group = if_else(interlayer.x != interlayer.y, 1, 0),
    ccode1a = if_else(ccode1 < ccode2, ccode1, ccode2),
    ccode2a = if_else(ccode1 > ccode2, ccode1, ccode2),
    ccode1 = ccode1a,
    ccode2, ccode2a
  ) %>% distinct(
    ccode1, ccode2, year, .keep_all = T
  )


full_regression_no_imputes <- fixest::feglm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio | ccode1 + ccode2 + dyad + year, data = dyad_data, family = "binomial")
full_lpm_no_imputes <- fixest::feols(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio | ccode1 + ccode2 + dyad + year, data = dyad_data)

full_regression_no_imputes_mid <- fixest::feglm(mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio | ccode1 + ccode2 + dyad + year, data = dyad_data, family = "binomial")
full_lpm_no_imputes_mid <- fixest::feols(mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio | ccode1 + ccode2 + dyad + year, data = dyad_data)


summary(full_regression_no_imputes_mid)
summary(full_lpm_no_imputes_mid)


save(full_regression_no_imputes, file = "model_output/full_regression_no_imputes.rda")
save(full_lpm_no_imputes, file = "model_output/full_lpm_no_imputes.rda")
save(full_regression_no_imputes_mid, file = "model_output/full_regression_no_imputes_mid.rda")
save(full_lpm_no_imputes_mid, file = "model_output/full_lpm_no_imputes_mid.rda")

gc()
logit_regression_fp_tau        <- glm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + taub + poly(year, 2), data = dyad_data, family = "binomial")
logit_regression_fp_kappaba    <- glm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + kappava  + poly(year, 2), data = dyad_data, family = "binomial")
lpm_regression_fp_tau          <- lm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + taub + poly(year, 2), data = dyad_data)
lpm_regression_fp_kappaba      <- lm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + kappava  + poly(year, 2), data = dyad_data)
logit_fe_regression_fp_tau     <- fixest::feglm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + taub | ccode1 + ccode2 + dyad + year, data = dyad_data, family = "binomial")
lpm_fe_fp_tau                  <- fixest::feols(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + taub | ccode1 + ccode2 + dyad + year, data = dyad_data)
logit_fe_regression_fp_kappaba <- fixest::feglm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + kappava | ccode1 + ccode2 + dyad + year, data = dyad_data, family = "binomial")
lpm_fe_fp_kappaba              <- fixest::feols(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + kappava | ccode1 + ccode2 + dyad + year, data = dyad_data)


logit_regression_fp_pi        <- glm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + piva + poly(year, 2), data = dyad_data, family = "binomial")
lpm_regression_fp_pi          <- lm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + piva + poly(year, 2), data = dyad_data)
logit_fe_fp_pi                <- fixest::feglm(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + piva | ccode1 + ccode2 + dyad + year, data = dyad_data, family = "binomial")
lpm_fe_fp_pi                  <- fixest::feols(fatal_mid ~ out_group + democ_joint + autoc_joint + major_power_joint + cinc_ratio + piva | ccode1 + ccode2 + dyad + year, data = dyad_data)


gc()
save(logit_regression_fp_tau, file = "model_output/logit_regression_fp_tau.rda")
save(logit_regression_fp_kappaba, file = "model_output/logit_regression_fp_kappaba.rda")
save(logit_regression_fp_pi, file = "model_output/logit_regression_fp_pi.rda")
save(lpm_regression_fp_tau, file = "model_output/lpm_regression_fp_tau.rda")
save(lpm_regression_fp_kappaba, file = "model_output/lpm_regression_fp_kappaba.rda")
save(lpm_regression_fp_pi, file = "model_output/lpm_regression_fp_pi.rda")
save(logit_fe_regression_fp_tau, file = "model_output/logit_fe_regression_fp_tau.rda")
save(logit_fe_regression_fp_kappaba, file = "model_output/logit_fe_regression_fp_kappaba.rda")
save(logit_fe_fp_pi, file = "model_output/logit_fe_fp_pi.rda")
save(lpm_fe_fp_tau, file = "model_output/lpm_fe_fp_tau.rda")
save(lpm_fe_fp_pi, file = "model_output/lpm_fe_fp_pi.rda")
save(lpm_fe_fp_kappaba, file = "model_output/lpm_fe_fp_kappaba.rda")


